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Abstract 

A two-dimensional lattice Boltzmann model with 19 discrete velocities for compressible fluids 
is proposed. The fifth-order Weighted Essentially Non-Oscillatory (5th-WENO) finite difference 
scheme is employed to calculate the convection term of the lattice Boltzmann equation. The validity 
of the model is verified by comparing simulation results of the Sod shock tube with its corresponding 
analytical solutions. The velocity and density gradient effects on the Kelvin-Helmholtz instability 
(KHI) are investigated using the proposed model. Sharp density contours are obtained in our 
simulations. It is found that, the linear growth rate 7 for the KHI decreases with increasing the 
width of velocity transition layer D v but increases with increasing the width of density transition 
layer D p . After the initial transient period and before the vortex has been well formed, the linear 
growth rates, 7„ and 7 P , vary with D v and D p approximately in the following way, In 7^ = a — bD v 
and 7 P = c + eln D p (D p < D p ), where a, b, c and e are fitting parameters, and D p is the effective 
interaction width of density transition layer. When D p > D p the linear growth rate 7^ does not 
vary significantly any more. One can use the hybrid effects of velocity and density transition layers 
to stabilize the KHI. Our numerical simulation results are in general agreement with the analytical 
results [L. F. Wang, et al, Phys. Plasma 17, 042103 (2010)]. 

PACS numbers: 47.11.-j, 47.20.4t, 05.20.Dd 

Keywords: lattice Boltzmann method; Kelvin-Helmholtz instability; velocity gradient effect; density gra- 
dient effect 
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I. INTRODUCTION 



As a mesoscopic approach and a bridge between the molecular dynamics method at 
the microscopic level and the conventional numerical method at the macroscopic level, the 
lattice Boltzmann (LB) method Q has been successfully applied to various fields during 
ihepast two decades, ranging from the multiphase system |2H4|, magnetohydrodynamics 
5-7], reaction-diffusion system jsljic|. compressible fluid dynamics 11 -18], simulations of 
linear and nonlinear partial differential equations 19] , etc. Its popularity is mainly owed 



to its kinetic nature 



20], which makes the physics at mesoscopic scale can be incorporated 



easily. As a result, this approach contains more physical connotation than Navier-Stokes 
or Euler equations based on the continuum hypothesis. Its popularity is also owed to its 
linear convection term, easy implementation of boundary conditions, simple programming 
and high efficiency in parallel computing, etc. 

The Kelvin- Helmholtz instability (KHI) occurs when two fluids flow with different tangen- 



tial velocities 



2l|. Under the condition of the KHI, small perturbations a 



between two fluids undergo linear 



22 



on g th e interface 



27H3l|. and may 



21426] and nonlinear growth stages 
evolve into turbulent mixing through nonlinear interactions. The KHI has attracted much 
attention because of its significance in both fundamental research and engineering applica- 



tions 



321 ] . On the one hand, in fundamental research, KHI is of great importance in the 



34 



and the interaction of the solar 



fields of turbulent mixing [33] , supernovae dynamics 

wind with the earth's magnetosphere {37], etc. On the other hand, in engineering applica- 



38 



40- 



and 



42j]. In the 



tions, the KHI plays a key role in small-scale mixing of Rayleigh- Taylor (K 
Richtmyer-Meshkov (RM) instabilities in inertial confinement fusion (ICF) 
final regime of RT and RM instabilities, KHI is initiated due to the shear velocity difference 
at the spike tips [43J] , and, therefore, the appearance of the KHI aggravates the development 
of final nonlinearity of RT instability or RM instability, and quickens the process of fluid 
flock mixing round the interface. 

In the present work, we propose a two-dimensional LB mode to simulate the KHI. The 
model consists of 19 discrete velocities in six directions and allows to recover the compress- 
ible Euler equations in the continuum limit. Actually, this phenomenon has been studied 



extensively by many researchers using experimental measurements^theoretical ana 



more recently by numerical simulations during the past decades |2lll25|, 



26 



31 



ysis and 



33]. Those 
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results indicate that the evolution of KHI depends on the density ratio, viscosities, sur- 
face tension, compressibility and others. Although the basic behavior of the KHI has been 
studied extensively, to the best of our knowledge, the use of the LB method to investigate 
the evolution of this phenomenon is still very limited. In this paper, with the LB method 
we focus on the velocity and the density gradient effects on KHI. The rest of the paper 
is organized as follows. In Section II, using the Chapman-Enskog analysis, we show that 
the current model can recover the Euler equations in the continuum limit. The numerical 
scheme and the validation of the model will be performed in Section III. In Section IV, we 
show numerical simulation results on the KHI for various widths of velocity and density 
transition layers. Both the two kinds of gradients effects are investigated carefully. Finally, 
in Section V conclusions and discussions are drawn. 



II. FORMULATION OF THE LB MODEL 



The model described here is inspired by the previous work of Watari and Tsutahara 



131 ] . which is based on a multispeed approach and where the truncated local equilibrium 



distribution uses global coefficients. 



A. Designing of the discrete velocity model 

In the present study, we will propose a finite difference LB (FDLB) model for compressible 
flows. Within this scheme, the physical symmetry of hydrodynamic systems can be more 
conveniently recovered by the use of discrete velocity model with high spatial isotropy. For 
completeness, let us start with the standard LB equation, an exact lattice-based dynamical 
equation expressed in terms of particle streaming and the Bhatnagar-Gross-Krook (BGK) 
collision operator 

fl(r + v;AM + At) - f' t (r,t) = -^(/; - if'), (1) 

where f- is the distribution function along the i direction, is its corresponding equilib- 
rium distribution function, r' is the spatial coordinate, is the particle velocity in the i 
direction, At' is the time step, and r is relaxation time determining the viscosity of fluid. 
For the standard LB method, all the particle velocities are restricted to those that exactly 
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link the lattice sites in unit time. In other words, there is a banding of discretizations of 
space and time, which is an inconvenience in the use of the standard LB model to study the 
thermohydrodynamic problems. 

The standard LB equation can be approximately expressed as the discrete velocity Boltz- 
mann equation 

Q t > ^ x ki q^i — r >\Jki Jki )■• \ z ) 

where v ki represents the discrete velocity model, subscript k indicates the kth group of the 
particle velocities whose speed is v' k . Eq.@ may be written in a dimensionless form by 
choosing three independent reference variables, the characteristic flow length scale L , the 
reference speed c s and the reference density p . The characteristic time to can be expressed 
as to = Lq/c s . With these reference variables, we have the following dimensionless variables 

* = ~, r = -, v H = — , /„ = — , t = 1 —. (3) 

c r c s Po to 

Consequently, we have the following dimensionless lattice Boltzmann equation 

To formulate a FDLB model, the next step is to chose an appropriate dimensionless 
discrete velocity model v fci . For a discrete velocity model described by 

v = 0, v ki = v k [cos{^-), sin(^-)], (5) 

the nth rank tensor is defined as 

n 
i=l 

where ati, «2,...,a n indicate either the x or y component. The tensor is isotropic if it is 
invariant under the coordinate rotation and the reflection. As for being isotropic, the odd 
rank tensors should be zero and the even rank tensors should be the sum of all possible 
products of Kronecker delta. In this study, the discrete velocity model with iV = 6 is used 
to construct a FDLB model. For this discrete velocity model it is easy find that the odd 
rank tensors are zero, and the even rank tensors generally have the following forms 

6 

2j v kiaVkip = 3vl5 a f3, (7) 
i=l 



VkiaVki/iVk^Vkix — ^ v k^a 



1=1 



(8) 



where 



Therefore, this discrete velocity model is isotropic up to, at least, its fifth rank tensor. 
B. Recovering of the Euler equations 



The macroscopic density p, momentum pu and temperature T are calculated as the 
moments of the local equilibrium distribution function 



/ 



(o) 

ki 



ki 



1 

1 , 



\ 



/ 



/ 



V 



p 

pu 

pT 



\ 



J 



(10) 



\ 2i V ki - U) 

Besides the moments in Eq. ffTO]) . Chapman- Enskog analysis indicates the following addi- 
tional ones are needed to satisfy in order to recover the hydrodynamic equations at Euler 
level 

^ VkiaVkipfif = P5 a p + pu a u p , (11) 



ki 



\ V li V kiafki = U *( 2P + \pU 2 )- 



ki 

To perform the Chap man- Enskog expansion on the two sides of Eq.( 
expansions 

Jki — Jki t Jki "r t Jki ' 



d d 2 d 

d d 
e 



Of 



dr r 



dr lc 



(12) 

we introduce 

(13) 
(14) 

(15) 



where e 1 is the Knudsen number. The introducing of e is equivalent to stipulating that 
the gas is dominated by large collision frequency. The fj^ is the zeroth order, fj£\ -^j- and 
o^- are the first order, and are the second order terms of the Knudsen number e. 
Substituting Eqs.( ll3p - (fl5"l) into Eq.(J3]) and comparing the coefficients of the same order of e 



give 



dt dr la 



T J ki ' 



(16) 



2 . Wki | d fki | d (Al).. A_ 1 f (2) / 17 n 

g /fc° ) | d fki | g/jS 1} | 5 rtU- 1 ),, \- 1 f(?) ns^ 



Summing Eq.lfl6l) over k,i gives 



Summing Eq. ffT8]) over k,i gives 



Using Eq.^HD and Eq.(fT5l) gives 



W + T-^ = - (19) 
dt dr la 



ih> 0,(j > 2). (20) 



| + ^(P«.) = 0. (2D 

It is clear that the continuity equation can be derived at any order of e. Performing the 
operator J2 v kia to the two sides of Eq. ffl6|) gives 



hi 



d d 

Tr-(pu a ) + ~ — [p(T5 a p + U a Up)} = 0. (22) 



Performing the operator Vkia to the two sides of Eq. ( Tl7l) gives 

(pu a ) = 0. (23) 



hi 

d_ 

Using Eq. (ll4p and Eq.( TT5|) gives the moment equation at the Euler level 



d d 

T^(pMa) + -^r[p(T5 a/3 + u a up)} = 0. (24) 



Performing operator to both sides of Eq. ffl6|) gives 

hi 



Performing operator Yl v ki/^ ^° both sides of Eq. fllT]) gives 

ki 

£«r + = (26) 
Using Eq.dHJ) and Eq.( TT5|) gives the energy equation at the Euler level 
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C. Formulation of the discrete equilibrium distribution function 



We now formulat the discrete equilibrium distribution function based on the discrete 
velocity model with N = Q. The requirement Eg. (1121) contains up to the third order of 
the flow velocity u. It is reasonable to expand the local equilibrium distribution fj^ in 
polynomial of the flow velocity up to the third order 

2 2 

= P*k[{i - -^) H j^^ 1 ~ H 

6T3 J ' { ' 



where 



'' exp[-J;] (29) 



2ttT " l 2T j 

is a function of temperature T and particle velocity v The truncated equilibrium distribu- 
tion function contains the third rank tensor of the particle velocity and the requirement 
Eg. ( ITT]) contains the second rank tensor. Thus, the discrete velocity model should have 
an isotropy up to its fifth rank tensor. So the discrete velocity model with N = 6 is an 
appropriate choice. 

To numerically calculate the equilibrium distribution function, one needs first calculate 
the global factor Fk. It should be noted that Fk can not be calculated directly from its 
definition Eq. (f29|) . It should take values in such a way that satisfies Eqs. f|T0|) -f fT2|) . 

To satisfy the first equation in Eq.( !T0|) . we require 

ki 

FkVkieVki V u £ u v = Tu 2 . (31) 

ki 

To satisfy the second equation in Eq. lflOl . we require 

ia^kie^e = Tu a , (32) 

ki 

FkVkiaVkieVkivVkH-UeUnUt = 3T 2 u 2 u a . (33) 

ki 

To satisfy the third equation in Eq.f JTOj) . we require 

E F 4 = T ' (34) 

ki 



vl 



^2 F k ^-v ki£ v kiv u £ u v = 2TV. (35) 

ki 

To satisfy Eq.( lTT|) . we require 

^ F k v kia v ki p = T5 a p, (36) 

ki 

^ FkVteaVkipVkisVkirtUeUT, = T 2 (u 2 5 al3 + 2u a u /3 ). (37) 
To satisfy Eg. ( 112]) . we require 

22,F k -tv kia v kie u e = 2T 2 u a , (38) 



hi 

i 

2 



^ Fk-^Vki a VkieVkir,Vk^U e U v U^ = 9T 3 u 2 u a . (39) 

ki 

If further consider the isotropic properties of the discrete velocity model, the above ten 
requirements reduce to four ones. Requirement Eq. (l30l) gives 



S>* = L (40) 



hi 



Requirements Eqs.flSI},©,®,® give 



E^4- ( 41 ) 

k 6 
Requirements Eqs. (133]) , (135]) , (137) , (|38|) give 

£^=|r 2 . ( 42 ) 

Requirement Eq. (l39|) gives 

J2F k vt=8T 3 . (43) 

To satisfy the above four requirements, four particle velocities v k are necessary. We choose 
a zero speed Vq = 0, and other three nonzero ones. Eqs.( l40|) - (|43|) are easily solved to the 
following solution 

Fk = 24T3 - 4(v 2 k+1 + v 2 +2 )T 2 + vJ ± A ± £ 
3v 2 (v 2 - v 2 k+1 )(v 2 - v 2 k+2 ) 
F = 1-6(F 1 + F 2 + F 3 ), (45) 

with 

.k + l, if k + I < 3; 
k + l = { ~ 1 = 1,2. (46) 

k + Z - 3, if k + I > 3. 
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The FDLB scheme breaks the binding of discretizations of space and time and makes 
the particle speeds more flexible. As far as a simulation being stably conducted, the specific 
values of v\, t>2 and V3 do not affect the accuracy of the results. The flexibility can be used to 
obtain the temperature range as wide as possible. In our simulations we set V\ — 1, V2 — 2 
and v 3 = 3. 

Moreover, it should be noted that the physical process described by the LB model has an 
intrinsic viscosity which is linearly related to the relaxation time r. The Euler equations or 
the Navier-Stokes equations are only approximations of the LB modeling in the continuum 
limit. In our case, we ensure only that the LB model recovers the correct Euler equations 
when the Knudsen number is zero. 



III. SPATIOTEMPORAL DISCRETIZATIONS AND VERIFICATION OF THE 
MODEL 



A. Time and space discretizations 

In this section, we will confirm validity of the model by conducting numerical simula- 
tions. Here the time derivative is calculated using the forward Euler FD scheme. Spatial 
derivatives in convection term Vki-jrfki are calculated using fifth order Weighted Essentially 



Non-Oscillatory (5th-WENO) FD scheme [44]. 



The WENO scheme is an improvement of the Essentially Non-Oscillatory (ENO) scheme, 
which changes the method of choosing smooth stencil with logical judgment into weighted 
average of all stencils, thereby improving the accuracy, the computational efficiency, and the 
adaptability for compressible flows containing shock waves, contact discontinuities, etc. 

The basic idea of the ENO scheme is to use the "smoothest" stencil among several candi- 
dates to approximate the numerical fluxes to a high order accuracy, and at the same time, to 
avoid spurious oscillations near shocks. The ENO scheme is uniformly high order accurate, 
robust, and there is essentially no oscillation near shock waves |45.]. However, this scheme 
has several drawbacks that can be improved by WENO scheme. For example, compared to 
WENO scheme, the ENO scheme is less accurate in smooth regions, and the numerical flux 
is not as smooth as that of the WENO scheme. For WENO scheme, instead of approximat- 
ing the numerical flux using only one of the candidate stencils, it uses all candidate stencils 
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through a convex combination. The contribution of each stencil to the final approximation 
of the numerical flux is determined by a weight, which can be defined in such a way that 
in smooth regions it approaches certain optimal weights to achieve a higher order of ac- 
curacy, while in regions near discontinuities, the stencils which contain the discontinuities 
are assigned a nearly zero weight [44J. As a result, the accuracy can be improved to the 



optimal order (in Ref. 46], a rth order ENO scheme has been improved to a (r + l)th 



WENO scheme, and in Ref. 44j the (r + l)th WENO scheme has been further improved to 
a (2r — l)th WENO scheme) in smooth regions and the essentially non-oscillatory property 
near discontinuities can be maintained. 

Now we spell out the details of the 5th- WENO scheme. To be specific, the discrete 
evolution equation of distribution function is 

rn+l _ rn d(. V kiafkij) . 1, n WOK*. 

JH,I — Jki,I a, ^ ~\Jki,I Jki,I J 1 ^ 1 

(If Q T 

= fu,i - 4r (W» " 'V-i/a) ^ - \Uki - (47) 

where the superscripts n, n + 1 indicate the consecutive two iteration steps, J is the index 
of lattice node, At is the time step, and ^7+1/2 is the numerical flux and can be defined as 
a convex combination 

hi,I+l/2 = Vih iiI+1 / 2 + ^2^7+1/2 + ^3^/+l/2- (48) 

Under the condition Vkia > 0, interpolation functions h q I+1 ^ 2 (q = 1, 2, 3) are given by 

1 7 11 

2 15 1 

K,I+l/2 — + g-^i,/ + (50) 

3 ' 6 ' + 6 
where F i}1 = v kia f kiJ . 

The weighting factors u q reflect the smooth degree of the template. They can be defined 
as follows 

oj = ^ q Co = Sq (52) 

q cD 1 + a} 2 + a} 3 ' q (10-6 + a,) 2 ' 1 1 

where 5\ = 1/10, 62 = 3/5 and £3 = 3/10 are the optimal weights. The small value 10~ 6 is 

added to the denominator to avoid dividing by zero. The coefficients a q in Eq. (l52l are the 
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fr|l+l/2 = oFiJ + ~ ^F i>I+2 , (51) 



smoothness indicators, and can be obtained by 

13 1 

ox = ^(Fi,J-2 - 2F i)/ _ 1 + Fij) 2 + - A {Fij-2 - AFij-! + 3Fu) 2 , (53) 

°2 = - 2F hI + F lJ+1 ) 2 + JCFi.j.! - P t , m ) 2 , (54) 

13 1 

a 3 = — (ify - 2F iJ+1 + F lJ+2 ) 2 + -(3F itI - AF hI+1 + F lJ+2 ) 2 . (55) 

Similarly, under the condition Vkia < 0, a mirror image procedure of the procedure from 
Eqs. (l48l to (|55l) can be carried out. 



B. Numerical test: one-dimensional Sod shock tube 

The validity of the formulated LB model is verified through a one-dimensional Riemann 

n 

problem, the Sod shock tube [47|. This is a classical test in the research of compressible 
flows, which contains the shock wave, the rarefaction wave and the contact discontinuity. 
For the problem considered, the initial condition is described by 

(p, it, v, T)\ L = (1.0, 0.0, 0.0, 1.0), x < 0; 

(56) 

(p, u, v, T)\ R = (0.125, 0.0, 0.0, 0.8), x > 0. 

Subscript "L" and U R" indicate macroscopic variables at the left and right sides of the 
discontinuity. The size of grid is Ax = Ay = 10~ 3 , time step and relaxation time are 
At = r = 10- 5 . 

The periodic boundary conditions are taken in the y direction. In the x direction, we set 

fki,-2,t = fki,-l,t — fki,0,t — /fei,l,t=o> (57) 

where —2, —1 and are the indexes of ghost nodes on the left side. Such a boundary 
condition means that the system at the boundaries keeps at their corresponding equilibrium 



states and the macroscopic quantities on the boundary nodes keep their initial values |17l.ll8|. 



(p, u, v, P)i=-2,t = (p, u, v, P)i=-i,t = (p, u, v, P)i=o, t = (p, u, v, P)t= m=0 . (58) 

On the right boundary, we can do in a similar way. The boundary conditions in the re- 
direction are also referred to as the supersonic inflow boundary conditions • Eq. (|57l) and 
Eq. (!58l) may be referred to as the microscopic and the macroscopic boundary conditions, 
respectively, which are consistent with each other. 
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FIG. 1: (Color online) Comparisons between LB results and exact ones for the one-dimensional 
Sod problem, where t = 0.2. 



It is worth noting that when the external environment is in equilibrium the boundary 
conditions for the FDLB model are essentially the same as those in traditional computa- 
tional fluid dynamics (CFD). When the external environment is not in equilibrium, the 
non-equilibrium part f^f 1 can be obtained from the inner lattice nodes via the extrapola- 
tion method. That nonequilibrium effects can be incorporated from boundary conditions is 
a merit of LB method over the traditional CFD. 

Figure 1 shows the computed density, pressure, velocity, and temperature profiles at 
t = 0.2, where the circles are for simulation results and solid lines are for analytical so- 
lutions. The two sets of results have a satisfying agreement. Figure 2 enlarges the part 
containing the shock wave for a closing view. The first-order upwind (lst-upwind) scheme, 
the second order upwind (2nd-upwind) scheme, the second order total variation diminishing 
(2nd-TVD) scheme 48[, and the 5th-WENO scheme are adopted in spatial discretization 
for comparisons. One can see that the lst-upwind scheme has a strong "smoothing effect" 
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0.385 0.390 0.395 0.400 0.385 0.390 0.395 0.400 
X X 

FIG. 2: (Color online) Local details of the profiles near the shock wave for Sod shock problem 
at t = 0.2. The lst-upwind scheme, the 2nd-upwind scheme, the 2nd-TVD scheme, and the 5th- 
WENO scheme are adopted in spatial discretization for comparisons. 

and gives a very smeared solution with excessive numerical dissipation. Compared to the 
former, the 2nd-upwind scheme has a weak "smoothing effect", but it introduces the unphys- 
ical oscillations at discontinuities. The 2nd-TVD scheme effectively eliminates the spurious 
numerical oscillations at discontinuities, while the numerical dissipation may need to be 
reduced further. However, with the 5th-WENO scheme, not only the spurious numerical 
oscillations at discontinuities are effectively refrained but also the numerical dissipation is 
severely curtailed. With an increase in the order of the schemes, fewer nodes are needed to 
capture the shock waves. The widths of the shock waves are spread over only about three 
to four grid cells with the 5th-WENO scheme, which shows that the LB model with the 
WENO scheme has a high ability to capture shocks for this test. 
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IV. EFFECTS OF VELOCITY AND DENSITY GRADIENTS ON KHI 



At the initial linear increasing stage of KHI, the amplitude r/ of perturbation evolves 
according to the following relation, rj = ^o e7 * ; where 7 is the growth coefficient and is 
dependent on the gradient of tangential velocity and gradient of density around the interface 



42| | . In other words, 7 is dependent on the width of velocity transition layer D v and width 
of density transition layer D p . In this section, we discuss separately the KHI in three cases, 
(i) D v is variable and D p is fixed, (ii) D p is variable and D v is fixed, (iii) both D p and D v 
are variable. The increasing rate 7 for case (i), (ii) and (iii) are referred to as j v , 7 P and 
7b, respectively We numerically obtain j v , j p and jr via fitting the curves of hiE x \ mgx (t) 
versus the time t, where E x \ max (t) is the maximum of E x (x, y, t) in the whole computational 
domain, E x (x, y, t) = p(x, y, t)u 2 (x, y, t)/2 is the perturbed kinetic energy at the position (x, 
y) at each time step t. 

A. Linear growth rate and velocity gradient effect 

The initial configurations in our simulation are described by 

p(*) = ^-^tanh(_l), (59) 

«M = ^-a tanh( JL), (60) 

Pl = Pr = P, (61) 

where D p and D v are the widths of density and velocity transition layers. The veloc- 
ity(density) is discontinuous at x = when D v = 0(D p = 0). Pl = 5.0(p^ = 2.0) is the 
density away from the interface of the left(right) fluid, vl = 0.5(ur = —0.5) is the velocity 
away from the interface in y-direction of the left(right) fluid and Pl(Pr)— 2.5 is the pressure 
in the left (right) side. The Mach numbers of the left side and the right side of the simulation 
domain are 0.5 and 0.32, respectively. Hence the flow speed vq is subsonic and the system 
can be approximately thought of as "incompressible". The whole calculation domain is a 
rectangle with length 0.6 and height 0.2, which is divided into 600 x 200 uniform meshes. A 
simple velocity perturbation in the x-direction is introduced to trigger the KH rollup and it 
is in the following form 

u = u sin(ky) exp(-kx), (62) 
15 




FIG. 3: (Color online) Density evolutions of KHI simulated using the LB model, where D v = 4 
and D p = 8,t = 0.1 in (a), t = 0.3 in (b), t = 0.5 in (c), and t = 0.7 in (d). 

where Uq = 0.02 is the amplitude of the perturbation. Here, k is the wave number of 
the initial perturbation, and is set to be 107T. The time step is At = 10~ 5 . Extensive 
studies indicate that viscosity will damp the evolution of the KHI. However, for the KHI in 
some cases, for example, in ICF, effects of the viscosity are generally negligible. Therefore, 
throughout the simulations, r is set to be 10~ 5 to reduce the physical viscosity. Boundary 
conditions for the simulations of KHI are as below. Periodic boundary conditions are used in 
the ^/-direction. In the horizontal direction, we adopt the outflow (zero gradient) boundary 



conditions, which are widely used by other authors to simulate the KHI [49rolj]. For the 
outflow boundary conditions, the zeroth-order extrapolation is used [17]]. According to this 
boundary conditions, at the left side, we set 

(p, U, V, P)i=-2,t = (p, u, v, P)i=-i, t = (p, u, v, P)i= ,t = {P, u, v, P)i=i,t- (63) 
and at the right side 

(p, U, V, P) I=Nx+3 ,t = (p, U, V, P)l= Nx +2,t = (P, U, V, P) J=Na+ i lt = (p, U, V, P)l=N x ,t, (64) 

where I = — 2, — 1, and N x + 1, N x + 2, N x + 3 are the indexes of left and the right ghost 
nodes, respectively. 

Figure 3 shows the temporal evolution of the density field for D v = 4 and D p = 8 at 
four different times. It is clear that at t = 0.3 the interface is wiggling due to the initial 
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FIG. 4: (Color online) Averaged density profiles along the x-axis at t = 0.1 (solid line), t = 0.3 
(dash dot line), t = 0.5 (dash dot dot line), and t = 0.7 (short dash line) for D p = 8 and D v = 4. 

perturbation and the velocity shear. After the initial linear growth stage, there is a nicely 
rolled vortex developing around the initial interface. A larger vortex is observed in the 
snapshot at t — 0.7, and a mixing layer is expected to be formed around the initial interface. 
The interface is continuous and smooth, which indicates the LB model has a good capturing 
ability of interface deformation. 

To quantitatively describe the characteristics of the vortex or the mixing layer, in Fig.4 
we plot the averaged density profile ~p{x) against the x-axis at t = 0.1, 0.3, 0.5 and 0.7. The 
averaged density in the mixing layer is defined as 



The averaged density profiles vary from being smooth to being irregular. The thickness of 
the mixing layer and the amplitude of the density oscillation increase with time. The zig- 
zags in the profiles indicate the transfer of fluids from the dense to the rarefactive regions 
and the irregularity in the mixing layer. 

To consider the velocity gradient effect, the simulations with various D v have been per- 
formed and the density maps for D v = 4, 8, 12, and 16 with five contour lines at t — 0.6 
are plotted in Fig. 5. We see that the width of the velocity transition layer can significantly 
affect the evolution of KHI. For fixed density gradient and velocity difference, the larger the 
value of D v , i.e., the wider the initial transition zone, the weaker the KHI, and the later 




(65) 
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(c) (d) 

FIG. 5: (Color online) Vortices in the mixing layer as a function of D v at t = 0.6, where D v = 4 in 
(a), D v = 8 in (b), D v = 12 in (c), and D v = 16 in (d). The density transition layer D p is fixed to 
be 8. 

the vortex appears. In Fig. 5(a) and Fig. 5(b), the appearance of large vortices demonstrates 
that the evolution is embarking on the nonlinear stage. While in Fig. 5(c) and Fig. 5(d), the 
width of the mixing layers is rapidly decreased by the increase of D v and the evolution is in 
the weakly nonlinear stage. Figures 5(a)-(d) show that a wider velocity transition zone has 
a better stabilization effect on KHI. 

Taking the logarithm of the perturbed peak kinetic energy E x \ ma/X of each time step, 
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FIG. 6: (Color online) Time evolution of the perturbed peak kinetic energy E x \ mauX along the x-axis 
in ln-linear scale for various widths of velocity transition layer. The dash-dotted lines represent 
the linear fits to the initial linear growth regimes. 
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FIG. 7: (Color online) Linear growth rate as a function of the width D v of velocity transition layer. 
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-|54[ (see Fig. 6). Therefore, 7 



the exponential growth manifests itself as a linear slope 
can be obtained from the slope of a linear function fitted to the initial growth stage. We 
note in this respect that E x oc u 2 oc (e 7 *) 2 grows at twice the rate of the KHI, namely 
k = 27, where k is the slope of the linear fit. The peak kinetic energy E x \ ma _ x can represent 
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the interacting strength of two different fluids. It is clear from Fig. 6 that the logarithm 
of the perturbed kinetic energy first grows linearly in time and then arrives at a nonlinear 
saturation procedure at late times. 

For a fixed D v and D p , E x \ max increases with time t during the linear growth stage. 
However, at the same time, the larger the value of D v , the smaller the perturbed peak 
kinetic energy E x \ max _. This means when the transition layer becomes wider, the evolution 
of the density field becomes slower. Moreover, we find that the dependence of j v on D v can 
be fitted by a logarithmic function with the form 

lnj v = a — bD v , (66) 

with a = 2.70 and b = 0.07 as displayed in Fig. 7. The numerical simulation result is in general 
agreement with the analytical results (see Eq.(18) and Fig. 3 in recent work of Wang, et al. 



26(|). In the classical case, the linear growth rate is 7 C = ky/pip 2 (vi — v 2 )/(pi + p 2 ) oc Av, 
where Av is the shear velocity difference. A wider transition layer decreases the local or the 
effective shear velocity difference Av, which results in a smaller linear growth rate, a smaller 
saturation energy, and a longer linear growth time tu n . 



B. Linear growth rate and density gradient effect 

In this subsection, the density gradient effect is investigated in a similar way. The initial 
conditions are described as, (pl, vl, Pl) = (5.0, 0.5, 1.5) and (pr, vr, Pr) = (1.25, —0.5, 
1.5). Parameters are set to be dx = dy = 0.002, At = 10~ 5 . Figure 8 shows time evolution 
of the logarithm of the peak kinetic energy E x \ m3X along the x-axis with various widths 
of density transition layers. It is seen in Fig. 8, for fixed width D v of velocity transition 
layer and fixed density difference, the linear growth rate first increases with the width D p of 
density transition layer. But when D p > 6, the linear growth rate does not vary significantly 
any more (see Fig.8 for details), which indicates the effective interaction width of D p is less 
than that of D v . Moreover, we find when D p < 6, the dependence of j p on D p can be fitted 
by the following equation 

7 P = c + e In D p , (67) 

with c = 5.28 and e = 0.62 as shown in the legend of Fig. 9. The numerical simulation results 
are also in general agreement with the analytical results (see Fig. 2 in previous work by Wang, 

20 




FIG. 8: (Color online) Time evolution of the logarithm of the peak kinetic energy E x 
the x-axis for various widths of density transition layers. 
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FIG. 9: (Color online) Linear growth rate as a function of the width D p of density transition layer. 

et al. 25], Eq.(18) and Fig.2 in recent work of Wang, et al. [26]). In the classical case, the 
square of the linear growth rate is 7^ = k 2 pip2(v\ — V2) 2 /(pi + P2) 2 oc (1 — v4 2 )At> 2 , where 
A = (pi — P2)/ (pi + P2) is the Atwood number. A wider density transition zone reduces the 
Atwood number around the interface. Then in the process of exchanging momentum in the 
direction normal to the interface, the perturbation can obtain more energy from the shear 
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FIG. 10: (Color online) The linear growth rate versus the width of density transition layer for 
R = 0.5, 1, 2, and 5. The initial density, shear velocity and pressure of the two fluids are {pi, vl, 
P L ) = (5.0, 0.5, 1.5) and (p R , v R , P R ) = (1.25, -0.5, 1.5). 

kinetic energy than in cases with sharper interfaces. Therefore, a wider density transition 
layer increases the linear growth rate of the KHI. 

C. Hybrid effects of velocity and density gradients 

In practical systems, at the interface of two fluids with a tangential velocity difference, 
both the velocity and the density gradients exist. The wider velocity transition layer de- 
creases the linear growth rate of the KHI, but the wider density transition layer strengthens 
it. Consequently, there is a competition between these two gradient effects. In this section, 
we investigate how these two effects compete with each other. For convenience, we introduce 
a coefficient R = D p /D v to analyze the combined effect. The linear growth rate at various 
velocity and density gradients for R = 0.5, 1, 2, and 5 are illustrated in Fig. 10. 

As shown in Fig. 10, on the whole, the hybrid effect of velocity and density transition 
layers reduces the linear growth rate 7^. Only at small D p and when R > 1, the hybrid 
effect of velocity and density transition layers makes larger the linear growth rate. This 
indicates that the effective interaction width of the velocity transition layer is wider 
than that of density transition layer Dp. 
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For a fixed value of R (> 1), when D p is small, the transfer efficiency of kinetic energy 
between neighboring layers increases with decreasing the density gradient, so the linear 
growth rate increases with D p and reaches a peak. When D p is greater than the peak value, 
the stabilizing effect of velocity transition layer becomes more significant and stronger than 
the destabilizing effect of the density transition layer, which leads to the decreasing of linear 
growth rate with D p , as shown in Fig. 10 of R = 1, 2, and 5 cases. However, when R < 1, 
the linear growth rate decreases monotonically with D p , no matter it is small or large. 

V. CONCLUSIONS AND DISCUSSIONS 

In this paper, a two-dimensional LB model with 19 discrete velocities in six directions is 
proposed. The model allows to recover the compressible Euler equations in the continuum 
limit. With the introduction of the 5th-WENO FD scheme, the unphysical oscillations at 
discontinuities are effectively suppressed and the numerical dissipation is severely curtailed. 
The validity of the model is verified by its application to the Sod shock tube, and excellent 
agreement between the simulation results and the analytical solutions can be found. 

Using the proposed LB model, the velocity and density gradient effects on the KHI 
have been studied. The averaged density profiles, and the peak kinetic energy are used to 
quantitatively describe the evolution of KHI. It is found that evolution of the KHI is damped 
with increasing the width of velocity transition layer but is strengthened with increasing the 
width of density transition layer. After the initial transient period and before the vortex has 
been well formed, the linear growth rates, 7^ and j p , vary with D v and D p approximately 
in the following way, ln7„ = a — bD v and 7 P = c + e\nD p (D p < D^), where a, b, c and e are 
fitting parameters and D p is the effective interaction width of density transition layer. When 
D p > the growth rate 7 P goes to nearly a constant and the destabilizing effect of density 
gradient on KHI is fixed. These can be understood as follows. As noted above, in the classical 
case, the square of the linear growth rate is 7^ = k 2 pip 2 (vi — v 2 ) 2 / (pi + P2) 2 oc (1 — A 2 )Av 2 . 
The wider velocity transition layer reduces the the local velocity difference Av and the wider 
density transition zone reduces the local Atwood number A. Therefore, the linear growth 
rate can be decreased by D v but made larger by D p . In practical system, both the two 
transition layers exist and there is a competition between these two gradient effects. One 
can use the hybrid effects of density and velocity transition layers to stabilize the KHI. By 
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incorporating an appropriate equation of state, or equivalently, a free energy functional, or 
an external force, the present model may be used to simulate the liquid-vapor transition and 
the surface tension effects on KHI. 
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